home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / bessel_olver.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  32.0 KB  |  983 lines
  1. /* specfunc/bessel_olver.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_airy.h>
  26.  
  27. #include "error.h"
  28.  
  29. #include "bessel.h"
  30. #include "bessel_olver.h"
  31.  
  32. #include "chebyshev.h"
  33. #include "cheb_eval.c"
  34.  
  35. /* fit for f(x) = zofmzeta((x+1)/2), 0 <= mzeta <= 1 */
  36. static double zofmzeta_a_data[20] = {
  37.   2.9332563730829348990,
  38.   0.4896518224847036624,
  39.   0.0228637617355380860,
  40.  -0.0001715731377284693,
  41.  -0.0000105927538148751,
  42.   1.0595602530419e-6,
  43.  -4.68016051691e-8,
  44.   5.8310020e-12,
  45.   1.766537581e-10,
  46.  -1.45034640e-11,
  47.   4.357772e-13,
  48.   4.60971e-14,
  49.  -2.57571e-14,
  50.   2.26468e-14,
  51.  -2.22053e-14,
  52.   2.08593e-14,
  53.  -1.84454e-14,
  54.   1.50150e-14,
  55.  -1.06506e-14,
  56.   5.5375e-15
  57. };
  58. static cheb_series zofmzeta_a_cs = {
  59.   zofmzeta_a_data,
  60.   19,
  61.   -1,1,
  62.   8
  63. };
  64.  
  65.  
  66. /* fit for f(x) = zofmzeta((9x+11)/2), 1 <= mzeta <= 10 */
  67. static double zofmzeta_b_data[30] = {
  68.   22.40725276466303489,
  69.   10.39808258825165581,
  70.   1.092050144486018425,
  71.  -0.071111274777921604,
  72.   0.008990125336059704,
  73.  -0.001201950338088875,
  74.   0.000106686807968315,
  75.   0.000017406491576830,
  76.  -0.000014946669657805,
  77.   6.189984487752e-6,
  78.  -2.049466715178e-6,
  79.   5.87189458020e-7,
  80.  -1.46077514157e-7,
  81.   2.9803936132e-8,
  82.  -3.817692108e-9,
  83.  -4.66980416e-10,
  84.   5.83860334e-10,
  85.  -2.78825299e-10,
  86.   1.01682688e-10,
  87.  -3.1209928e-11,
  88.   8.111122e-12,
  89.  -1.663986e-12,
  90.   1.81364e-13,
  91.   5.3414e-14,
  92.  -4.7234e-14,
  93.   2.1689e-14,
  94.  -7.815e-15,
  95.   2.371e-15,
  96.  -6.04e-16,
  97.   1.20e-16
  98. };
  99. static cheb_series zofmzeta_b_cs = {
  100.   zofmzeta_b_data,
  101.   29,
  102.   -1,1,
  103.   15
  104. };
  105.  
  106.  
  107. /* fit for f(x) = zofmzeta(mz(x))/mz(x)^(3/2),
  108.  * mz(x) = (2/(x+1))^(2/3) 10
  109.  * 10 <= mzeta <= Inf
  110.  */
  111. static double zofmzeta_c_data[11] = {
  112.   1.3824761227122911500,
  113.   0.0244856101686774245,
  114.  -0.0000842866496282540,
  115.   1.4656076569771e-6,
  116.  -3.14874099476e-8,
  117.   7.561134833e-10,
  118.  -1.94531643e-11,
  119.   5.245878e-13,
  120.  -1.46380e-14,
  121.   4.192e-16,
  122.  -1.23e-17
  123. };
  124. static cheb_series zofmzeta_c_cs = {
  125.   zofmzeta_c_data,
  126.   10,
  127.   -1,1,
  128.   6
  129. };
  130.  
  131.  
  132. /* Invert [Abramowitz+Stegun, 9.3.39].
  133.  * Assumes minus_zeta >= 0.
  134.  */
  135. double
  136. gsl_sf_bessel_Olver_zofmzeta(double minus_zeta)
  137. {
  138.   if(minus_zeta < 1.0) {
  139.     const double x = 2.0*minus_zeta - 1.0;
  140.     gsl_sf_result c;
  141.     cheb_eval_e(&zofmzeta_a_cs, x, &c);
  142.     return c.val;
  143.   }
  144.   else if(minus_zeta < 10.0) {
  145.     const double x = (2.0*minus_zeta - 11.0)/9.0;
  146.     gsl_sf_result c;
  147.     cheb_eval_e(&zofmzeta_b_cs, x, &c);
  148.     return c.val;
  149.   }
  150.   else {
  151.     const double TEN_32 = 31.62277660168379332; /* 10^(3/2) */
  152.     const double p = pow(minus_zeta, 3.0/2.0);
  153.     const double x = 2.0*TEN_32/p - 1.0;
  154.     gsl_sf_result c;
  155.     cheb_eval_e(&zofmzeta_c_cs, x, &c);
  156.     return c.val * p;
  157.   }
  158. }
  159.  
  160.  
  161. /* Chebyshev fit for f(x) = z(x)^6 A_3(z(x)),  z(x) = 22/(10(x+1)) */
  162. static double A3_gt1_data[31] = {
  163.   -0.123783199829515294670493131190,
  164.    0.104636462534700704670877382304,
  165.   -0.067500816575851826744877535903,
  166.    0.035563362418888483652711005520,
  167.   -0.0160738524035979408472979609051,
  168.    0.0064497878252851092073278056238,
  169.   -0.00235408261133449663958121821593,
  170.    0.00079545702851302155411892534965,
  171.   -0.00025214920745855079895784825637,
  172.    0.00007574004596069392921153301833,
  173.   -0.00002172917966339623434407978263,
  174.    5.9914810727868915476543145465e-06,
  175.   -1.5958781571808992162953719817e-06,
  176.    4.1232986512903717525448312012e-07,
  177.   -1.0369725993417659101913919101e-07,
  178.    2.5457982304266541145999235022e-08,
  179.   -6.1161715053791743082427422443e-09,
  180.    1.4409346199138658887871461320e-09,
  181.   -3.3350445956255561668232014995e-10,
  182.    7.5950686572918996453336138108e-11,
  183.   -1.7042296334409430377389900278e-11,
  184.    3.7723525020626230919721640081e-12,
  185.   -8.2460237635733980528416501227e-13,
  186.    1.7816961527997797696251868875e-13,
  187.   -3.8084101506541792942694560802e-14,
  188.    8.0593669930916099079755351563e-15,
  189.   -1.6896565961641739017452636964e-15,
  190.    3.5115651805888443184822853595e-16,
  191.   -7.2384771938569255638904297651e-17,
  192.    1.4806598977677176106283840244e-17,
  193.   -3.0069285750787303634897997963e-18
  194. };
  195. static cheb_series A3_gt1_cs = {
  196.   A3_gt1_data,
  197.   30,
  198.   -1,1,
  199.   17
  200. };
  201.  
  202. /* chebyshev expansion for f(x) = z(x)^8 A_4(z(x)), z(x) = 12/(5(x+1)) */
  203. static double A4_gt1_data[30] = {
  204.   1.15309329391198493586724229008,
  205.  -1.01812701728669338904729927846,
  206.   0.71964022270555684403652781941,
  207.  -0.42359963977172689685150061355,
  208.   0.215024488759339557817435404261,
  209.  -0.096751915348145944032096342479,
  210.   0.039413982058824310099856035361,
  211.  -0.014775225692561697963781115014,
  212.   0.005162114514159370516947823271,
  213.  -0.00169783446445524322560925166335,
  214.   0.00052995667873006847211519193478,
  215.  -0.00015802027574996477115667974856,
  216.   0.000045254366680989687988902825193,
  217.  -0.000012503722965474638015488600967,
  218.   3.3457656998119148699124716204e-06,
  219.  -8.6981575241150758412492331833e-07,
  220.   2.2030895484325645640823940625e-07,
  221.  -5.4493369492600677068285936533e-08,
  222.   1.3190457281724829107139385556e-08,
  223.  -3.1301560183377379158951191769e-09,
  224.   7.2937802527123344842593076131e-10,
  225.  -1.6712080137945140407348940109e-10,
  226.   3.7700053248213600430503521194e-11,
  227.  -8.3824538848817227637828899571e-12,
  228.   1.8388741910049766865274037194e-12,
  229.  -3.9835919980753778560117573063e-13,
  230.   8.5288827136546615604290389711e-14,
  231.  -1.8060227869114416998653266836e-14,
  232.   3.7849342199690728470461022877e-15,
  233.  -7.8552867468122209577151823365e-16
  234. };
  235. static cheb_series A4_gt1_cs = {
  236.   A4_gt1_data,
  237.   17, /* 29, */
  238.   -1, 1,
  239.   17
  240. };
  241.  
  242. /* Chebyshev fit for f(x) = z(x)^3 B_2(z(x)), z(x) = 12/(5(x+1)) */
  243. static double B2_gt1_data[40] = {
  244.   0.00118587147272683864479328868589,
  245.   0.00034820459990648274622193981840,
  246.  -0.00030411304425639768103075864567,
  247.   0.00002812066284012343531484682886,
  248.   0.00004493525295901613184489898748,
  249.  -0.00003037629997093072196779489677,
  250.   0.00001125979647123875721949743970,
  251.  -2.4832533969517775991951008218e-06,
  252.  -9.9003813640537799587086928278e-08,
  253.   4.9259859656183110299492296029e-07,
  254.  -3.7644120964426705960749504975e-07,
  255.   2.2887828521334625189639122509e-07,
  256.  -1.3202687370822203731489855050e-07,
  257.   7.7019669092537400811434860763e-08,
  258.  -4.6589706973010511603890144294e-08,
  259.   2.9396476233013923711978522963e-08,
  260.  -1.9293230611988282919101954538e-08,
  261.   1.3099107013728717842406906896e-08,
  262.  -9.1509111940885962831104149355e-09,
  263.   6.5483472971925614347299375295e-09,
  264.  -4.7831253582139967461241674569e-09,
  265.   3.5562625457426178152760148639e-09,
  266.  -2.6853389444008414186916562103e-09,
  267.   2.0554738667134200145781857289e-09,
  268.  -1.5923172019517426277886522758e-09,
  269.   1.2465923213464381457319481498e-09,
  270.  -9.8494846881180588507969988989e-10,
  271.   7.8438674499372126663957464312e-10,
  272.  -6.2877567918342950225937136855e-10,
  273.   5.0662318868755257959686944117e-10,
  274.  -4.0962270881243451160378710952e-10,
  275.   3.3168684677374908553161911299e-10,
  276.  -2.6829406619847450633596163305e-10,
  277.   2.1603988122184568375561077873e-10,
  278.  -1.7232373309560278402012124481e-10,
  279.   1.3512709089611470626617830434e-10,
  280.  -1.0285354732538663013167579792e-10,
  281.   7.4211345443901713467637018423e-11,
  282.  -4.8124980266864320351456993068e-11,
  283.   2.3666534694476306077416831958e-11
  284. };
  285. static cheb_series B2_gt1_cs = {
  286.   B2_gt1_data,
  287.   39,
  288.   -1, 1,
  289.   30
  290. };
  291.  
  292.  
  293. /* Chebyshev fit for f(x) = z(x)^6 B_3(z(x)), z(x) = 12/(5(x+1)) */
  294. static double B3_gt1_data[30] = {
  295.  -0.0102445379362695740863663926486,
  296.   0.0036618484329295342954730801917,
  297.   0.0026154252498599303282569321117,
  298.  -0.0036187389410353156728771706336,
  299.   0.0021878564157692275944613452462,
  300.  -0.0008219952303590803584426516821,
  301.   0.0001281773889155631494321316520,
  302.   0.0001000944653368032985720548637,
  303.  -0.0001288293344663774273453147788,
  304.   0.00010136264202696513867821487205,
  305.  -0.00007000275849659556221916572733,
  306.   0.00004694886396757430431607955146,
  307.  -0.00003190003869717837686356945696,
  308.   0.00002231453668447775219665947479,
  309.  -0.00001611102197712439539300336438,
  310.   0.00001196634424990735214466633513,
  311.  -9.0986920398931223804111374679e-06,
  312.   7.0492613694235423068926562567e-06,
  313.  -5.5425216624642184684300615394e-06,
  314.   4.4071884714230296614449244106e-06,
  315.  -3.5328595506791663127928952625e-06,
  316.   2.84594975572077091520522824686e-06,
  317.  -2.29592697828824392391071619788e-06,
  318.   1.84714740375289956396370322228e-06,
  319.  -1.47383331248116454652025598620e-06,
  320.   1.15687781098593231076084710267e-06,
  321.  -8.8174688524627071175315084910e-07,
  322.   6.3705856964426840441434605593e-07,
  323.  -4.1358791499961929237755474814e-07,
  324.   2.0354151158738819867477996807e-07
  325. };
  326. static cheb_series B3_gt1_cs = {
  327.   B3_gt1_data,
  328.   29,
  329.   -1, 1,
  330.   29
  331. };
  332.  
  333.  
  334. /* Chebyshev fit for f(x) = z(x) B_2(z(x)), z(x) = 2(x+1)/5 */
  335. static double B2_lt1_data[40] = {
  336.   0.00073681565841337130021924199490,
  337.   0.00033803599647571227535304316937,
  338.  -0.00008251723219239754024210552679,
  339.  -0.00003390879948656432545900779710,
  340.   0.00001961398056848881816694014889,
  341.  -2.35593745904151401624656805567e-06,
  342.  -1.79055017080406086541563835433e-06,
  343.   1.33129571185610681090725934031e-06,
  344.  -5.38879444715436544130673956170e-07,
  345.   1.49603056041381416881299945557e-07,
  346.  -1.83377228267274327911131293091e-08,
  347.  -1.33191430762944336526965187651e-08,
  348.   1.60642096463700438411396889489e-08,
  349.  -1.28932576330421806740136816643e-08,
  350.   9.6169275086179165484403221944e-09,
  351.  -7.1818502280703532276832887290e-09,
  352.   5.4744009217215145730697754561e-09,
  353.  -4.2680446690508456935030086136e-09,
  354.   3.3941665009266174865683284781e-09,
  355.  -2.7440714072221673882163135170e-09,
  356.   2.2488361522108255229193038962e-09,
  357.  -1.8638240716608748862087923337e-09,
  358.   1.5592350940805373500866440401e-09,
  359.  -1.3145743937732330609242633070e-09,
  360.   1.1153716777215047842790244968e-09,
  361.  -9.5117576805266622854647303110e-10,
  362.   8.1428799553234876296804561100e-10,
  363.  -6.9893770813548773664326279169e-10,
  364.   6.0073113636087448745018831981e-10,
  365.  -5.1627434258513453901420776514e-10,
  366.   4.4290993195074905891788459756e-10,
  367.  -3.7852978599966867611179315200e-10,
  368.   3.2143959338863177145307610452e-10,
  369.  -2.7025926680620777594992221143e-10,
  370.   2.2384857772457918539228234321e-10,
  371.  -1.8125071664276678046551271701e-10,
  372.   1.4164870008713668767293008546e-10,
  373.  -1.0433101857132782485813325981e-10,
  374.   6.8663910168392483929411418190e-11,
  375.  -3.4068313177952244040559740439e-11
  376. };
  377. static cheb_series B2_lt1_cs = {
  378.   B2_lt1_data,
  379.   39,
  380.   -1, 1,
  381.   39
  382. };
  383.  
  384.  
  385. /* Chebyshev fit for f(x) = B_3(2(x+1)/5) */
  386. static double B3_lt1_data[40] = {
  387.  -0.00137160820526992057354001614451,
  388.  -0.00025474937951101049982680561302,
  389.   0.00024762975547895881652073467771,
  390.   0.00005229657281480196749313930265,
  391.  -0.00007488354272621512385016593760,
  392.   0.00001416880012891046449980449746,
  393.   0.00001528986060172183690742576230,
  394.  -0.00001668672297078590514293325326,
  395.   0.00001061765189536459018739585094,
  396.  -5.8220577442406209989680801335e-06,
  397.   3.3322423743855900506302033234e-06,
  398.  -2.23292405803003860894449897815e-06,
  399.   1.74816651036678291794777245325e-06,
  400.  -1.49581306041395051804547535093e-06,
  401.   1.32759146107893129050610165582e-06,
  402.  -1.19376077392564467408373553343e-06,
  403.   1.07878303863211630544654040875e-06,
  404.  -9.7743335011819134006676476250e-07,
  405.   8.8729318903693324226127054792e-07,
  406.  -8.0671146292125665050876015280e-07,
  407.   7.3432860378667354971042255937e-07,
  408.  -6.6897926072697370325310483359e-07,
  409.   6.0966619703735610352576581485e-07,
  410.  -5.5554095284507959561958605420e-07,
  411.   5.0588335673197236002812826526e-07,
  412.  -4.6008146297767601862670079590e-07,
  413.   4.1761348515688145911438168306e-07,
  414.  -3.7803230006989446874174476515e-07,
  415.   3.4095248501364300041684648230e-07,
  416.  -3.0603959751354749520615015472e-07,
  417.   2.7300134179365690589640458993e-07,
  418.  -2.4158028250762304756044254231e-07,
  419.   2.1154781038298751985689113868e-07,
  420.  -1.8269911328756771201465223313e-07,
  421.   1.5484895085808513749026173074e-07,
  422.  -1.2782806851555809369226440495e-07,
  423.   1.0148011725394892565174207341e-07,
  424.  -7.5658969771439627809239950461e-08,
  425.   5.0226342286491286957075289622e-08,
  426.  -2.5049645660259882970547555831e-08
  427. };
  428. static cheb_series B3_lt1_cs = {
  429.   B3_lt1_data,
  430.   39,
  431.   -1, 1,
  432.   39
  433. };
  434.  
  435.  
  436. /* Chebyshev fit for f(x) = A_3(9(x+1)/20) */
  437. static double A3_lt1_data[40] = {
  438.   -0.00017982561472134418587634980117,
  439.   -0.00036558603837525275836608884064,
  440.   -0.00002819398055929628850294406363,
  441.    0.00016704539863875736769812786067,
  442.   -0.00007098969970347674307623044850,
  443.   -8.4470843942344237748899879940e-06,
  444.    0.0000273413090343147765148014327150,
  445.   -0.0000199073838489821681991178018081,
  446.    0.0000100004176278235088881096950105,
  447.   -3.9739852013143676487867902026e-06,
  448.    1.2265357766449574306882693267e-06,
  449.   -1.88755584306424047416914864854e-07,
  450.   -1.37482206060161206336523452036e-07,
  451.    2.10326379301853336795686477738e-07,
  452.   -2.05583778245412633433934301948e-07,
  453.    1.82377384812654863038691147988e-07,
  454.   -1.58130247846381041027699152436e-07,
  455.    1.36966982725588978654041029615e-07,
  456.   -1.19250280944620257443805710485e-07,
  457.    1.04477169029350256435316644493e-07,
  458.   -9.2064832489437534542041040184e-08,
  459.    8.1523798290458784610230199344e-08,
  460.   -7.2471794980050867512294061891e-08,
  461.    6.4614432955971132569968860233e-08,
  462.   -5.7724095125560946811081322985e-08,
  463.    5.1623107567436835158110947901e-08,
  464.   -4.6171250746798606260216486042e-08,
  465.    4.1256621998650164023254101585e-08,
  466.   -3.6788925543159819135102047082e-08,
  467.    3.2694499457951844422299750661e-08,
  468.   -2.89125899697964696586521743928e-08,
  469.    2.53925288725374047626589488217e-08,
  470.   -2.20915707933726481321465184207e-08,
  471.    1.89732166352720474944407102940e-08,
  472.   -1.60058977893259856012119939554e-08,
  473.    1.31619294542205876946742394494e-08,
  474.   -1.04166651771938038563454275883e-08,
  475.    7.7478015858156185064152078434e-09,
  476.   -5.1347942579352613057675111787e-09,
  477.    2.5583541594586723967261504321e-09
  478. };
  479. static cheb_series A3_lt1_cs = {
  480.   A3_lt1_data,
  481.   39,
  482.   -1, 1,
  483.   39
  484. };
  485.  
  486. /* chebyshev fit for f(x) = A_4(2(x+1)/5) */
  487. static double A4_lt1_data[30] = {
  488.   0.00009054703770051610946958226736,
  489.   0.00033066000498098017589672988293,
  490.   0.00019737453734363989127226073272,
  491.  -0.00015490809725932037720034762889,
  492.  -0.00004514948935538730085479280454,
  493.   0.00007976881782603940889444573924,
  494.  -0.00003314566154544740986264993251,
  495.  -1.88212148790135672249935711657e-06,
  496.   0.0000114788756505519986352882940648,
  497.  -9.2263039911196207101468331210e-06,
  498.   5.1401128250377780476084336340e-06,
  499.  -2.38418218951722002658891397905e-06,
  500.   1.00664292214481531598338960828e-06,
  501.  -4.23224678096490060264249970540e-07,
  502.   2.00132031535793489976535190025e-07,
  503.  -1.18689501178886741400633921047e-07,
  504.   8.7819524319114212999768013738e-08,
  505.  -7.3964150324206644900787216386e-08,
  506.   6.5780431507637165113885884236e-08,
  507.  -5.9651053193022652369837650411e-08,
  508.   5.4447762662767276209052293773e-08,
  509.  -4.9802057381568863702541294988e-08,
  510.   4.5571368194694340198117635845e-08,
  511.  -4.1682117173547642845382848197e-08,
  512.   3.8084701352766049815367147717e-08,
  513.  -3.4740302885185237434662649907e-08,
  514.   3.1616557064701510611273692060e-08,
  515.  -2.8685739487689556252374879267e-08,
  516.   2.5923752117132254429002796600e-08,
  517.  -2.3309428552190587304662883477e-08
  518. };
  519. static cheb_series A4_lt1_cs = {
  520.   A4_lt1_data,
  521.   29,
  522.   -1, 1,
  523.   29
  524. };
  525.  
  526.  
  527. static double olver_B0(double z, double abs_zeta)
  528. {
  529.   if(z < 0.98) {
  530.     const double t = 1.0/sqrt(1.0-z*z);
  531.     return -5.0/(48.0*abs_zeta*abs_zeta) + t*(-3.0 + 5.0*t*t)/(24.0*sqrt(abs_zeta));
  532.   }
  533.   else if(z < 1.02) {
  534.     const double a = 1.0-z;
  535.     const double c0 =  0.0179988721413553309252458658183;
  536.     const double c1 =  0.0111992982212877614645974276203;
  537.     const double c2 =  0.0059404069786014304317781160605;
  538.     const double c3 =  0.0028676724516390040844556450173;
  539.     const double c4 =  0.0012339189052567271708525111185;
  540.     const double c5 =  0.0004169250674535178764734660248;
  541.     const double c6 =  0.0000330173385085949806952777365;
  542.     const double c7 = -0.0001318076238578203009990106425;
  543.     const double c8 = -0.0001906870370050847239813945647;
  544.     return c0 + a*(c1 + a*(c2 + a*(c3 + a*(c4 + a*(c5 + a*(c6 + a*(c7 + a*c8)))))));
  545.   }
  546.   else {
  547.     const double t = 1.0/(z*sqrt(1.0 - 1.0/(z*z)));
  548.     return -5.0/(48.0*abs_zeta*abs_zeta) + t*(3.0 + 5.0*t*t)/(24.0*sqrt(abs_zeta));
  549.   }
  550. }
  551.  
  552.  
  553. static double olver_B1(double z, double abs_zeta)
  554. {
  555.   if(z < 0.88) {
  556.     const double t   = 1.0/sqrt(1.0-z*z);
  557.     const double t2  = t*t;
  558.     const double rz  = sqrt(abs_zeta);
  559.     const double z32 = rz*rz*rz;
  560.     const double z92 = z32*z32*z32;
  561.     const double term1 = t*t*t * (30375.0 - 369603.0*t2 + 765765.0*t2*t2 - 425425.0*t2*t2*t2)/414720.0;
  562.     const double term2 = 85085.0/(663552.0*z92);
  563.     const double term3 = 385.0/110592.*t*(3.0-5.0*t2)/(abs_zeta*abs_zeta*abs_zeta);
  564.     const double term4 = 5.0/55296.0*t2*(81.0 - 462.0*t2 + 385.0*t2*t2)/z32;
  565.     return -(term1 + term2 + term3 + term4)/rz;
  566.   }
  567.   else if(z < 1.12) {
  568.     const double a = 1.0-z;
  569.     const double c0  = -0.00149282953213429172050073403334;
  570.     const double c1  = -0.00175640941909277865678308358128;
  571.     const double c2  = -0.00113346148874174912576929663517;
  572.     const double c3  = -0.00034691090981382974689396961817;
  573.     const double c4  =  0.00022752516104839243675693256916;
  574.     const double c5  =  0.00051764145724244846447294636552;
  575.     const double c6  =  0.00058906174858194233998714243010;
  576.     const double c7  =  0.00053485514521888073087240392846;
  577.     const double c8  =  0.00042891792986220150647633418796;
  578.     const double c9  =  0.00031639765900613633260381972850;
  579.     const double c10 =  0.00021908147678699592975840749194;
  580.     return c0+a*(c1+a*(c2+a*(c3+a*(c4+a*(c5+a*(c6+a*(c7+a*(c8+a*(c9+a*c10)))))))));
  581.   }
  582.   else {
  583.     const double t   = 1.0/(z*sqrt(1.0 - 1.0/(z*z)));
  584.     const double t2  = t*t;
  585.     const double rz  = sqrt(abs_zeta);
  586.     const double z32 = rz*rz*rz;
  587.     const double z92 = z32*z32*z32;
  588.     const double term1 = -t2*t * (30375.0 + 369603.0*t2 + 765765.0*t2*t2 + 425425.0*t2*t2*t2)/414720.0;
  589.     const double term2 = 85085.0/(663552.0*z92);
  590.     const double term3 = -385.0/110592.0*t*(3.0+5.0*t2)/(abs_zeta*abs_zeta*abs_zeta);
  591.     const double term4 = 5.0/55296.0*t2*(81.0 + 462.0*t2 + 385.0*t2*t2)/z32;
  592.     return (term1 + term2 + term3 + term4)/rz;
  593.   }
  594. }
  595.  
  596.  
  597. static double olver_B2(double z, double abs_zeta)
  598. {
  599.   if(z < 0.8) {
  600.     const double x = 5.0*z/2.0 - 1.0;
  601.     gsl_sf_result c;
  602.     cheb_eval_e(&B2_lt1_cs, x, &c);
  603.     return  c.val / z;
  604.   }
  605.   else if(z <= 1.2) {
  606.     const double a = 1.0-z;
  607.     const double c0 = 0.00055221307672129279005986982501;
  608.     const double c1 = 0.00089586516310476929281129228969;
  609.     const double c2 = 0.00067015003441569770883539158863;
  610.     const double c3 = 0.00010166263361949045682945811828;
  611.     const double c4 = -0.00044086345133806887291336488582;
  612.     const double c5 = -0.00073963081508788743392883072523;
  613.     const double c6 = -0.00076745494377839561259903887331;
  614.     const double c7 = -0.00060829038106040362291568012663;
  615.     const double c8 = -0.00037128707528893496121336168683;
  616.     const double c9 = -0.00014116325105702609866850307176;
  617.     return c0+a*(c1+a*(c2+a*(c3+a*(c4+a*(c5+a*(c6+a*(c7+a*(c8+a*c9))))))));
  618.   }
  619.   else {
  620.     const double zi = 1.0/z;
  621.     const double x  = 12.0/5.0 * zi - 1.0;
  622.     gsl_sf_result c;
  623.     cheb_eval_e(&B2_gt1_cs, x, &c);
  624.     return c.val * zi*zi*zi;
  625.   }
  626. }
  627.  
  628.  
  629. static double olver_B3(double z, double abs_zeta)
  630. {
  631.   if(z < 0.8) {
  632.     const double x = 5.0*z/2.0 - 1.0;
  633.     gsl_sf_result c;
  634.     cheb_eval_e(&B3_lt1_cs, x, &c);
  635.     return c.val;
  636.   }
  637.   else if(z < 1.2) {
  638.     const double a = 1.0-z;
  639.     const double c0 = -0.00047461779655995980754441833105;
  640.     const double c1 = -0.00095572913429464297452176811898;
  641.     const double c2 = -0.00080369634512082892655558133973;
  642.     const double c3 = -0.00000727921669154784138080600339;
  643.     const double c4 =  0.00093162500331581345235746518994;
  644.     const double c5 =  0.00149848796913751497227188612403;
  645.     const double c6 =  0.00148406039675949727870390426462;
  646.     return c0 + a*(c1 + a*(c2 + a*(c3 + a*(c4 + a*(c5 + a*c6)))));
  647.   }
  648.   else {
  649.     const double x   = 12.0/(5.0*z) - 1.0;
  650.     const double zi2 = 1.0/(z*z);
  651.     gsl_sf_result c;
  652.     cheb_eval_e(&B3_gt1_cs, x, &c);
  653.     return  c.val * zi2*zi2*zi2;
  654.   }
  655. }
  656.  
  657.  
  658. static double olver_A1(double z, double abs_zeta, double * err)
  659. {
  660.   if(z < 0.98) {
  661.     double t = 1.0/sqrt(1.0-z*z);
  662.     double rz = sqrt(abs_zeta);
  663.     double t2 = t*t;
  664.     double term1 =  t2*(81.0 - 462.0*t2 + 385.0*t2*t2)/1152.0;
  665.     double term2 = -455.0/(4608.0*abs_zeta*abs_zeta*abs_zeta);
  666.     double term3 =  7.0*t*(-3.0 + 5.0*t2)/(1152.0*rz*rz*rz);
  667.     *err = 2.0 * GSL_DBL_EPSILON * (fabs(term1) + fabs(term2) + fabs(term3));
  668.     return term1 + term2 + term3;
  669.   }
  670.   else if(z < 1.02) {
  671.     const double a = 1.0-z;
  672.     const double c0 = -0.00444444444444444444444444444444;
  673.     const double c1 = -0.00184415584415584415584415584416;
  674.     const double c2 =  0.00056812076812076812076812076812;
  675.     const double c3 =  0.00168137865661675185484709294233;
  676.     const double c4 =  0.00186744042139000122193399504324;
  677.     const double c5 =  0.00161330105833747826430066790326;
  678.     const double c6 =  0.00123177312220625816558607537838;
  679.     const double c7 =  0.00087334711007377573881689318421;
  680.     const double c8 =  0.00059004942455353250141217015410;
  681.     const double sum = c0+a*(c1+a*(c2+a*(c3+a*(c4+a*(c5+a*(c6+a*(c7+a*c8)))))));
  682.     *err = 2.0 * GSL_DBL_EPSILON * fabs(sum);
  683.     return sum;
  684.   }
  685.   else {
  686.     const double t = 1.0/(z*sqrt(1.0 - 1.0/(z*z)));
  687.     const double rz = sqrt(abs_zeta);
  688.     const double t2 = t*t;
  689.     const double term1 = -t2*(81.0 + 462.0*t2 + 385.0*t2*t2)/1152.0;
  690.     const double term2 =  455.0/(4608.0*abs_zeta*abs_zeta*abs_zeta);
  691.     const double term3 = -7.0*t*(3.0 + 5.0*t2)/(1152.0*rz*rz*rz);
  692.     *err = 2.0 * GSL_DBL_EPSILON * (fabs(term1) + fabs(term2) + fabs(term3));
  693.     return term1 + term2 + term3;
  694.   }
  695. }
  696.  
  697.  
  698. static double olver_A2(double z, double abs_zeta)
  699. {
  700.   if(z < 0.88) {
  701.     double t  = 1.0/sqrt(1.0-z*z);
  702.     double t2 = t*t;
  703.     double t4 = t2*t2;
  704.     double t6 = t4*t2;
  705.     double t8 = t4*t4;
  706.     double rz = sqrt(abs_zeta);
  707.     double z3 = abs_zeta*abs_zeta*abs_zeta;
  708.     double z32 = rz*rz*rz;
  709.     double z92 = z3*z32;
  710.     double term1 = t4*(4465125.0 - 94121676.0*t2 + 349922430.0*t4 - 446185740.0*t6  + 185910725.0*t8)/39813120.0;
  711.     double term2 = -40415375.0/(127401984.0*z3*z3);
  712.     double term3 = -95095.0/15925248.0*t*(3.0-5.0*t2)/z92;
  713.     double term4 = -455.0/5308416.0 *t2*(81.0 - 462.0*t2 + 385.0*t4)/z3;
  714.     double term5 = -7.0/19906560.0*t*t2*(30375.0 - 369603.0*t2  + 765765.0*t4  - 425425.0*t6)/z32;
  715.     return term1 + term2 + term3 + term4 + term5;
  716.   }
  717.   else if(z < 1.12) {
  718.     double a = 1.0-z;
  719.     const double c0  =  0.000693735541354588973636592684210;
  720.     const double c1  =  0.000464483490365843307019777608010;
  721.     const double c2  = -0.000289036254605598132482570468291;
  722.     const double c3  = -0.000874764943953712638574497548110;
  723.     const double c4  = -0.001029716376139865629968584679350;
  724.     const double c5  = -0.000836857329713810600584714031650;
  725.     const double c6  = -0.000488910893527218954998270124540;
  726.     const double c7  = -0.000144236747940817220502256810151;
  727.     const double c8  =  0.000114363800986163478038576460325;
  728.     const double c9  =  0.000266806881492777536223944807117;
  729.     const double c10 = -0.011975517576151069627471048587000;
  730.     return c0+a*(c1+a*(c2+a*(c3+a*(c4+a*(c5+a*(c6+a*(c7+a*(c8+a*(c9+a*c10)))))))));
  731.   }
  732.   else {
  733.     const double t  = 1.0/(z*sqrt(1.0 - 1.0/(z*z)));
  734.     const double t2 = t*t;
  735.     const double t4 = t2*t2;
  736.     const double t6 = t4*t2;
  737.     const double t8 = t4*t4;
  738.     const double rz = sqrt(abs_zeta);
  739.     const double z3 = abs_zeta*abs_zeta*abs_zeta;
  740.     const double z32 = rz*rz*rz;
  741.     const double z92 = z3*z32;
  742.     const double term1 = t4*(4465125.0 + 94121676.0*t2 + 349922430.0*t4 + 446185740.0*t6  + 185910725.0*t8)/39813120.0;
  743.     const double term2 = -40415375.0/(127401984.0*z3*z3);
  744.     const double term3 =  95095.0/15925248.0*t*(3.0+5.0*t2)/z92;
  745.     const double term4 = -455.0/5308416.0 *t2*(81.0 + 462.0*t2 + 385.0*t4)/z3;
  746.     const double term5 =  7.0/19906560.0*t*t2*(30375.0 + 369603.0*t2  + 765765.0*t4  + 425425.0*t6)/z32;
  747.     return term1 + term2 + term3 + term4 + term5;
  748.   }
  749. }
  750.  
  751.  
  752. static double olver_A3(double z, double abs_zeta)
  753. {
  754.   if(z < 0.9) {
  755.     const double x = 20.0*z/9.0 - 1.0;
  756.     gsl_sf_result c;
  757.     cheb_eval_e(&A3_lt1_cs, x, &c);
  758.     return c.val;
  759.   }
  760.   else if(z < 1.1) {
  761.     double a = 1.0-z;
  762.     const double c0 = -0.000354211971457743840771125759200;
  763.     const double c1 = -0.000312322527890318832782774881353;
  764.     const double c2 =  0.000277947465383133980329617631915;
  765.     const double c3 =  0.000919803044747966977054155192400;
  766.     const double c4 =  0.001147600388275977640983696906320;
  767.     const double c5 =  0.000869239326123625742931772044544;
  768.     const double c6 =  0.000287392257282507334785281718027;
  769.     return c0 + a*(c1 + a*(c2 + a*(c3 + a*(c4 + a*(c5 + a*c6)))));
  770.   }
  771.   else {
  772.     const double x   = 11.0/(5.0*z) - 1.0;
  773.     const double zi2 = 1.0/(z*z);
  774.     gsl_sf_result c;
  775.     cheb_eval_e(&A3_gt1_cs, x, &c);
  776.     return  c.val * zi2*zi2*zi2;
  777.   }
  778. }
  779.  
  780.  
  781. static double olver_A4(double z, double abs_zeta)
  782. {
  783.   if(z < 0.8) {
  784.     const double x = 5.0*z/2.0 - 1.0;
  785.     gsl_sf_result c;
  786.     cheb_eval_e(&A4_lt1_cs, x, &c);
  787.     return c.val;
  788.   }
  789.   else if(z < 1.2) {
  790.     double a = 1.0-z;
  791.     const double c0 =  0.00037819419920177291402661228437;
  792.     const double c1 =  0.00040494390552363233477213857527;
  793.     const double c2 = -0.00045764735528936113047289344569;
  794.     const double c3 = -0.00165361044229650225813161341879;
  795.     const double c4 = -0.00217527517983360049717137015539;
  796.     const double c5 = -0.00152003287866490735107772795537;
  797.     return c0 + a*(c1 + a*(c2 + a*(c3 + a*(c4 + a*c5))));
  798.   }
  799.   else {
  800.     const double x   = 12.0/(5.0*z) - 1.0;
  801.     const double zi2 = 1.0/(z*z);
  802.     gsl_sf_result c;
  803.     cheb_eval_e(&A4_gt1_cs, x, &c);
  804.     return c.val * zi2*zi2*zi2*zi2;
  805.   }
  806. }
  807.  
  808. inline
  809. static double olver_Asum(double nu, double z, double abs_zeta, double * err)
  810. {
  811.   double nu2 = nu*nu;
  812.   double A1_err;
  813.   double A1 = olver_A1(z, abs_zeta, &A1_err);
  814.   double A2 = olver_A2(z, abs_zeta);
  815.   double A3 = olver_A3(z, abs_zeta);
  816.   double A4 = olver_A4(z, abs_zeta);
  817.   *err = A1_err/nu2 + GSL_DBL_EPSILON;
  818.   return 1.0 + A1/nu2 + A2/(nu2*nu2) + A3/(nu2*nu2*nu2) + A4/(nu2*nu2*nu2*nu2);
  819. }
  820.  
  821. inline
  822. static double olver_Bsum(double nu, double z, double abs_zeta)
  823. {
  824.   double nu2 = nu*nu;
  825.   double B0 = olver_B0(z, abs_zeta);
  826.   double B1 = olver_B1(z, abs_zeta);
  827.   double B2 = olver_B2(z, abs_zeta);
  828.   double B3 = olver_B3(z, abs_zeta);
  829.   return B0 + B1/nu2 + B2/(nu2*nu2) + B3/(nu2*nu2*nu2*nu2);
  830. }
  831.  
  832.  
  833. /* uniform asymptotic, nu -> Inf, [Abramowitz+Stegun, 9.3.35]
  834.  *
  835.  * error:
  836.  *    nu =  2: uniformly good to >  6D
  837.  *    nu =  5: uniformly good to >  8D
  838.  *    nu = 10: uniformly good to > 10D
  839.  *    nu = 20: uniformly good to > 13D
  840.  *
  841.  */
  842. int gsl_sf_bessel_Jnu_asymp_Olver_e(double nu, double x, gsl_sf_result * result)
  843. {
  844.   /* CHECK_POINTER(result) */
  845.  
  846.   if(x <= 0.0 || nu <= 0.0) {
  847.     DOMAIN_ERROR(result);
  848.   }  
  849.   else {
  850.     double zeta, abs_zeta;
  851.     double arg;
  852.     double pre;
  853.     double asum, bsum, asum_err;
  854.     gsl_sf_result ai;
  855.     gsl_sf_result aip;
  856.     double z = x/nu;
  857.     double crnu = pow(nu, 1.0/3.0);
  858.     double nu3  = nu*nu*nu;
  859.     double nu11 = nu3*nu3*nu3*nu*nu;
  860.     int stat_a, stat_ap;
  861.  
  862.     if(fabs(1.0-z) < 0.02) {
  863.       const double a = 1.0-z;
  864.       const double c0 = 1.25992104989487316476721060728;
  865.       const double c1 = 0.37797631496846194943016318218;
  866.       const double c2 = 0.230385563409348235843147082474;
  867.       const double c3 = 0.165909603649648694839821892031;
  868.       const double c4 = 0.12931387086451008907;
  869.       const double c5 = 0.10568046188858133991;
  870.       const double c6 = 0.08916997952268186978;
  871.       const double c7 = 0.07700014900618802456;
  872.       pre = c0 + a*(c1 + a*(c2 + a*(c3 + a*(c4 + a*(c5 + a*(c6 + a*c7))))));
  873.       zeta = a * pre;
  874.       pre  = sqrt(2.0*sqrt(pre/(1.0+z)));
  875.       abs_zeta = fabs(zeta);
  876.     }
  877.     else if(z < 1.0) {
  878.       double rt   = sqrt(1.0 - z*z);
  879.       abs_zeta = pow(1.5*(log((1.0+rt)/z) - rt), 2.0/3.0);
  880.       zeta = abs_zeta;
  881.       pre  = sqrt(2.0*sqrt(abs_zeta/(rt*rt)));
  882.     }
  883.     else {
  884.       /* z > 1 */
  885.       double rt = z * sqrt(1.0 - 1.0/(z*z));
  886.       abs_zeta = pow(1.5*(rt - acos(1.0/z)), 2.0/3.0);
  887.       zeta = -abs_zeta;
  888.       pre  = sqrt(2.0*sqrt(abs_zeta/(rt*rt)));
  889.     }
  890.  
  891.     asum = olver_Asum(nu, z, abs_zeta, &asum_err);
  892.     bsum = olver_Bsum(nu, z, abs_zeta);
  893.  
  894.     arg  = crnu*crnu * zeta;
  895.     stat_a  = gsl_sf_airy_Ai_e(arg, GSL_MODE_DEFAULT, &ai);
  896.     stat_ap = gsl_sf_airy_Ai_deriv_e(arg, GSL_MODE_DEFAULT, &aip);
  897.  
  898.     result->val  = pre * (ai.val*asum/crnu + aip.val*bsum/(nu*crnu*crnu));
  899.     result->err  = pre * (ai.err * fabs(asum/crnu));
  900.     result->err += pre * fabs(ai.val) * asum_err / crnu;
  901.     result->err += pre * fabs(ai.val * asum) / (crnu*nu11);
  902.     result->err += 8.0 * GSL_DBL_EPSILON * fabs(result->val);
  903.  
  904.     return GSL_ERROR_SELECT_2(stat_a, stat_ap);
  905.   }
  906. }
  907.  
  908.  
  909. /* uniform asymptotic, nu -> Inf,  [Abramowitz+Stegun, 9.3.36]
  910.  *
  911.  * error:
  912.  *    nu =  2: uniformly good to >  6D
  913.  *    nu =  5: uniformly good to >  8D
  914.  *    nu = 10: uniformly good to > 10D
  915.  *    nu = 20: uniformly good to > 13D
  916.  */
  917. int gsl_sf_bessel_Ynu_asymp_Olver_e(double nu, double x, gsl_sf_result * result)
  918. {
  919.   /* CHECK_POINTER(result) */
  920.  
  921.   if(x <= 0.0 || nu <= 0.0) {
  922.     DOMAIN_ERROR(result);
  923.   }  
  924.   else {
  925.     double zeta, abs_zeta;
  926.     double arg;
  927.     double pre;
  928.     double asum, bsum, asum_err;
  929.     gsl_sf_result bi;
  930.     gsl_sf_result bip;
  931.     double z = x/nu;
  932.     double crnu = pow(nu, 1.0/3.0);
  933.     double nu3  = nu*nu*nu;
  934.     double nu11 = nu3*nu3*nu3*nu*nu;
  935.     int stat_b, stat_d;
  936.  
  937.     if(fabs(1.0-z) < 0.02) {
  938.       const double a = 1.0-z;
  939.       const double c0 = 1.25992104989487316476721060728;
  940.       const double c1 = 0.37797631496846194943016318218;
  941.       const double c2 = 0.230385563409348235843147082474;
  942.       const double c3 = 0.165909603649648694839821892031;
  943.       const double c4 = 0.12931387086451008907;
  944.       const double c5 = 0.10568046188858133991;
  945.       const double c6 = 0.08916997952268186978;
  946.       const double c7 = 0.07700014900618802456;
  947.       pre = c0 + a*(c1 + a*(c2 + a*(c3 + a*(c4 + a*(c5 + a*(c6 + a*c7))))));
  948.       zeta = a * pre;
  949.       pre  = sqrt(2.0*sqrt(pre/(1.0+z)));
  950.       abs_zeta = fabs(zeta);
  951.     }
  952.     else if(z < 1.0) {
  953.       double rt   = sqrt(1.0 - z*z);
  954.       abs_zeta = pow(1.5*(log((1.0+rt)/z) - rt), 2.0/3.0);
  955.       zeta = abs_zeta;
  956.       pre  = sqrt(2.0*sqrt(abs_zeta/(rt*rt)));
  957.     }
  958.     else {
  959.       /* z > 1 */
  960.       double rt = z * sqrt(1.0 - 1.0/(z*z));
  961.       double ac = acos(1.0/z);
  962.       abs_zeta = pow(1.5*(rt - ac), 2.0/3.0);
  963.       zeta = -abs_zeta;
  964.       pre  = sqrt(2.0*sqrt(abs_zeta)/rt);
  965.     }
  966.  
  967.     asum = olver_Asum(nu, z, abs_zeta, &asum_err);
  968.     bsum = olver_Bsum(nu, z, abs_zeta);
  969.  
  970.     arg  = crnu*crnu * zeta;
  971.     stat_b = gsl_sf_airy_Bi_e(arg, GSL_MODE_DEFAULT, &bi);
  972.     stat_d = gsl_sf_airy_Bi_deriv_e(arg, GSL_MODE_DEFAULT, &bip);
  973.  
  974.     result->val  = -pre * (bi.val*asum/crnu + bip.val*bsum/(nu*crnu*crnu));
  975.     result->err  =  pre * (bi.err * fabs(asum/crnu));
  976.     result->err +=  pre * fabs(bi.val) * asum_err / crnu;
  977.     result->err +=  pre * fabs(bi.val*asum) / (crnu*nu11);
  978.     result->err +=  8.0 * GSL_DBL_EPSILON * fabs(result->val);
  979.  
  980.     return GSL_ERROR_SELECT_2(stat_b, stat_d);
  981.   }
  982. }
  983.